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In this paper, the flow in the diffuser section of the Icing Research Tunnel at 
the NASA Lewis Research Center is numerically investigated. To accomplish 
this, an existing computer code is utilized. The code, known as PARC3D, is 
based on the Beam-Warming algorithm applied to the strong conservation 
law form of the complete Navier-Stokes equations. The first portion of the 
paper consists of a brief description of the diffuser and its current flow 
characteristics. A brief discussion of the code work follows. Predicted velocity 
patterns are then compared with the measured values. 


INTRODUCTION 

The diffuser section of the Icing Research 
Tunnel (IRT) at the NASA Lewis Research Center 
has been used in the past to perform a limited 
amount of testing. However, measurements have 
shown that the flow quality is highly nonuniform 
and is unsuitable for many aerodynamic tests. This 
poor flow quality probably adversely affects the 
overall performance of the the IRT although the 
airflow in the test section appears to be well 
behaved. Because of this it was decided that there 
was a need to better understand the peculiar flow 
characteristics in the diffuser of the IRT. This 
understanding may suggest methods which may be 
used to improve the flow quality. If this objective is 
realized, the diffuser can be more widely used as a 
second, larger, but lower speed, test section. Also, 
the study may help to improve the overall 
performance of the tunnel. 

Diffusers are used to decelerate flows such 
that static pressure can be recovered and, in the 
case of wind tunnels, friction losses can be reduced. 
A number of studies have been performed to gather 
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information on performance and characteristics of 
diffusers. However, much of the fundamental data 
in the literature is from axisymmetric and/or two- 
dimensional diffuser studies. Little information 
exists for fully three-dimensional diffusers. The IRT 
diffuser, having four flat walls which all diverge at 
constant angles, is of this three-dimensional type. 

The study discussed herein is a 
continuation of the work reported in Ref. 1. That 
paper details the flow measurements which were 
made in the IRT to establish quantitatively some 
of the flow characteristics of the tunnel and to form 
a basis of comparison for computational studies of 
the flow. The initial computational results which 
precipitated the current numerical simulation are 
also discussed. 

Numerical simulation techniques are 
increasingly being used to model and study internal 
flows. Two such techniques have been used to 
model the flow in the diffuser of the IRT. The first 
technique, which has been incorporated into a 
computer code known as PEPSIG, uses a forward 
marching procedure to solve a parabolized form of 
the Navier-Stokes equations. This procedure has the 
capability of solving three-dimensional, turbulent, 
subsonic flow problems. The PEPSIG code has had 



success in the past in accurately modeling internal 
flows in inlets and ducts. However, it has some 
restrictions which limit its usefulness when applied 
to the IRT geometry. One limitation concerns its 
inability to model the rectangular shaped cross 
section of the IRT. Instead it used a superelliptical 
cross section which rounds each corner. This is 
similar to adding fillets to the cross section of the 
tunnel. Another limitation concerns its limited grid 
size, particularly in the axial direction. This 
restricts the amount of the test section and diffuser 
which could be accurately modelled. The flow 
predicted by the PEPSIG code was markedly 
different from that which was measured in the 
tunnel. The predicted velocity profiles in the 
diffuser were much flatter and more uniform than 
those measured. Also, the code computed a much 
smaller flow speed in the four corners of the tunnel 
than the measured values. PEPSIG predicted that 
some separation was occuring in the corners 
whereas the measurements strongly suggest that 
separation, if it does occur, would not occur in the 
corners, as will be discussed later. 

Due to the lack of agreement between the 
computed results from PEPSIG and the measured 
values, a second numerical simulation was 
employed. Application and discussion of this second 
computer program, known as the PARC3D code, is 
the subject of this paper. 


VENT DOOR 


IRT GEOMETRY AND MEASURED FLOW 
CHARACTERISTICS 

A diagram of the IRT is shown in Fig. 1. 
It is a dosed loop wind tunnel currently capable of 
operating at test section airspeeds of up to 300 
mph. A large refrigeration system allows the tunnel 
to operate at total temperatures as low as minus 
20*F. Air flows from a settling chamber through a 
nozzle with a contraction ratio of 14.13:1.0 into the 
test section which is 20 feet long and has a 
rectangular cross section that is 6 feet high by 9 
feet wide. The settling chamber and contraction 
section are also rectangular in shape. 

The portion of the IRT which forms the 
diffuser is denoted in Fig. 1. This portion is 81.5 
feet long and consists of four flat, straight walls 
each of which diverges from the centerline at an 
angle of 2.5*. The diffuser entrance is 6 feet high 
and 9 feet wide while the exit is 13.469 feet high 
and 16.469 feet wide. The expansion ratio of the 
diffuser is 4.11:1.0. 

As mentioned earlier, the flow in the test 
section and the diffuser entrance appears to be well 
behaved. Figure 2 shows measured vertical velocity 
profiles at the centerline of the diffuser entrance for 
test section speeds of 150 and 300 mph. Figure 3 
shows measured vertical velocity profiles in five 
equally spaced locations across the diffuser entrance 
for a test section speed of 200 mph. In each of these 
figures the vertical location of the local measure- 



Fig. 1 Icing Research Wind Tunnel = diffuser entrance, (£) —diffuser exit) 
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Fig. 2 Measured Profiles at Diffuser Entrance 

ment above the floor of the tunnel is plotted 
against the local value of the measurement which is 
nominalized by the tunnel reference value as 
measured by the facilty Pitot-static tube located in 
the test section. As may be seen, each of these 
profiles is relatively flat and uniform. 



■ .w V.U ■ U.o I .U U.o I .U I ,U 

NORTH N-CENTER CENTER S-CENTER SOUTH 
VELOCITY RATIO, V!oc /Vr*1 


Fig. 3 Measured Profiles Across Diffuser Entrance 

Flow in the diffuser exit is markedly 
different from that measured in the diffuser inlet. 
Figure 4 shows measured vertical velocity profiles 
at the centerline of the diffuser exit for test section 
airspeeds of 150 and 300 mph. As can be seen, 
these profiles are nonuniform. Figure 5 shows 
measured vertical velocity profiles at five equally 
spaced locations across the diffuser exit for a test 
section airspeed of 200 mph. Not only are these 
profiles nonuniform, but adjacent profiles are 
dissimilar. The north and south profiles in Fig. 5 
indicate that the flow is tending to migrate toward 
each of the four corners of the diffuser exit. 


Fig. 4 Measured Profiles at Diffuser Exit 
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Fig. 5 Measured Profiles Across Diffuser Exit 


The contour plot of measured local velocities at the 
diffuser exit at 150 mph is shown in Fig. 6 and also 
demonstrates this tendency in the flow. 

I«T DIFFUSES VELOCITY CONTOURS 



Fig. 6 Measured Contours at Diffuser Exit 
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It is interesting to also note that a 
similarity exists between the measured airflow 
velocity contour plot in the diffuser exit and a test 
section map of peak liquid water content from each 
spray bar in the system 5 . Figure 7 shows such a 
map of this peak LWC. This map was established 
by accreting ice on 2-inch diameter vertical bars 
mounted in the test section. Each of the eight 
horizontally mounted spray bars located in the 
settling chamber of the tunnel were individually 
operated. The horizontal lines numbered 1-8 in the 
map indicate the line of maximum ice accretion on 
each of the vertical bars in the test section. The 
numbered vertical lines in the map indicate the 
maximum ice accretion when a single vertical 
column of nozzles was operated. This similarity 
indicates the need for further aerodynamic 
investigation and suggests that these flow 
characteristics may not be limited to the diffuser. 



Fig. 7 Map of Peak LWC at Test Section 
THE PARC3D COMPUTER CODE 


An existing computer code 2 was used to 
simulate the flow in the IRT. This code is known 
as PARC3D. It is based on the complete Navier- 
Stokes equations written in strong conservation law 
form. The code uses the Beam and Warming 
approximate factorization algorithm 3 to solve a set 
of finite difference equations which were produced 
by central differencing the Navier-Stokes equations 
on a regular grid. The code calculates flow 
characteristics based on a specified boundary 
geometry and the corresponding flow conditions on 
these boundaries. A wide range of geometries and 
boundary conditions may be specified. 


The flow can be either inviscid or viscous. 
If viscous, it can be treated as either laminar or 
turbulent. The turbulence model used is the 
relatively simple algebraic model of Baldwin and 
Lomax 4 . 

It has been found that the code has had 
less success when the flow velocity is below a Mach 
number of 0.1. Moreover, the code does not usually 
give good results if the grid has high aspect ratio 
cells. Convergence of the PARC code can be 
adversely affected by tight grid packing. Use of the 
turbulence model also tends to slow convergence. 

CODE SET-UP FOR IRT SIMULATION 

The PARC3D code requires as input: a 
grid, initial and boundary conditions, and program 
execution controls. The initial condition, which 
may also be regarded as an initial solution, does 
not need to be an accurate one because only steady 
state results, which are obtained by an asymptotic 
time march, are of interest. Thus, only a rough, 
approximate solution is required to get the 
processing started. 

The primary boundary conditions used for 
the IRT flow simulation runs were: four no-slip, 
adiabatic walls and subsonic inlet and outlet flows. 
Uniform distributions of total pressure and 
temperature were specified on the inflow boundary. 
At the exit to the diffuser, a uniform static 
pressure, corresponding to a particular tunnel 
airspeed, was prescribed. The Reynolds number 
also had to be specified. For the test section 
airspeed of 150 mph, the Reynolds number was 
1.65xl0 6 . The program controls allow the user 
limited control over program execution such as: 
limiting the total number of iterations, maximum 
time step size, tabular output, and so forth. 

Several grids were created to reflect the 
IRT geometry. The initial grid developed for the 
IRT assumed a one-quarter plane symmetry in the 
cross section of the tunnel. However, because 
asymmetric effects had been observed in the tunnel 
flow, it was decided to use a full cross section 
m ode l of the tunnel. Such a cross section is shown 
in Fig. 8. The grid generation scheme developed for 
this work allowed for grid packing near the walls. 
This was desired in order to resolve the larger 
gradients in the flow that were anticipated in these 
regions. The grids generated had a cross section of 
29 points by 29 points and included 55 points in 
the axial direction. It is recognized that this grid is 
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a bit coarse, and may not have sufficient grid lines 
within the wall layers. However, it was deemed an 
adequate starting structure for this problem. 

Two three-dimensional configurations of 
the tunnel were used. The first was set up to 
include the aft one-half of the test section and the 
entire length of the diffuser as shown in Fig. 9. In 
order to study the effects of the contraction section 
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Fig. 8 Full Cross Section for Grid 



Fig. 9 Mesh for Half Test Section and Diffuser 

on diffuser flow, a second configuration was set up 
which included the settling chamber, the 
contraction section, the entire length of the test 
section, and the entire length of the diffuser. This 
grid is depicted in Fig. 10. 



Fig. 10 Mesh for Extended Grid 


RESULTS AND DISCUSSION 


The results from each particular 
combination of geometry and flow conditions will 
be presented through the use of both velocity 
contour plots and velocity profile plots at key 
locations along the length of the tunnel. The 
velocity profile plots provide a more detailed view 
of the predicted flow characteristics at a particular 
location while the velocity contour plots give an 
overall picture of the flow characteristics at the 
particular cross section of the tunnel. 

Even though the main leg of the IRT is 
geometrically symmetric in both the horizontal and 
vertical planes, it was felt that the numerical 
simulations should be run across the entire cross 
section rather than simply a one-quarter or one-half 
cross section. In this way, asymmetries due to wall 
obstacles or flow effects may be subsequently 
incorporated into the simulation. 

The first case was that of turbulent flow 
for a test section speed of 150 mph. The grid used 
includes the aft one-half of the test section and the 
entire diffuser as shown in Fig. 9. The grid 
dimensions were 29x29x55. Figure 11 shows the 
velocity contours at the first axial station of the 
grid which corresponds to a cross section midway 
down the test section of the tunnel. The contour 
levels are in terms of Mach numbers with the 
highest level corresponding to the centermost 
contour line. As may be seen, the contour lines are 
bunched relatively tightly together indicating a 
small boundary layer, and, as might be expected, 
are very regular and symmetric around the tunnel 
walls. The boundary layer is slightly thicker at 
both side walls than at the floor and ceiling. 
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Fig. 11 Predicted Velocity Contours, Test Section 
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Figure 12 shows predicted velocity contours 
at the entrance to the diffuser. Again the contours 
are regular and symmetric around the tunnel walls. 
The boundary layer has grown slightly toward the 
center of the tunnel and remains slightly thicker 
along the side walls. Figure 13 is the predicted 
vertical velocity profile at the tunnel centerline. 
The vectors are located at points on the grid and 
are proportional the the velocity at each point. 
These profiles also reflect the uniformity of the 
predicted flow at this location. 



Fig. 12 Predicted Vel. Contours, Diffuser Entrance 



Fig. 13 Predicted Vel. Profile, Diffuser Entrance 

Predicted velocity contours at the diffuser 
exit are shown in Fig. 14. They indicate a flow 
which is uniformly varying across the tunnel with a 
small amount of separation in or backflow occuring 
in each corner. Figure 15 shows the vertical velocity 
profile at the tunnel centerline. It also depicts the 
uniformly varying velocity predicted at this 
location. 

To investigate the effects of the contraction 
section on the flow in the diffuser, the grid was 
extended to include this part of the tunnel. This 
includes approximately 10 ft. of settling chamber, 
the 35 ft. of contraction section and the first 10 feet 
of the test section. The flow conditions were set to 
laminar flow for a test section airspeed of 150 mph. 



Fig. 14 Predicted Vel. Contours, Diffuser Exit 



Fig. 15 Predicted Vel. Profile, Diffuser Exit 

Figure 16 shows predicted velocity contours 
for these conditions at a station halfway down the 
test section. This corresponds to the initial station 
in the previous case. The contours are regular and 
symmetric around the walls of the tunnel. The 
boundary layer is slightly thicker here than in the 
previous case. 



Fig. 16 Predicted Vel. Contour, Test Section, 
Extended Grid, Laminar 


Figure 17 shows predicted velocity contours 
at the diffuser entrance. Again, the contours are 
regular and symmetric around the tunnel walls. 
Figure 18 presents a computed vertical velocity 
profile at the centerline of the diffuser entrance 
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Fig. 17 Predicted Vel. Contours, Diffuser Entrance, 
Extended Grid, Laminar 



Fig. 18 Predicted Vel. Profile, Diffuser Entrance, 
Extended Grid, Laminar 

which also indicates the well behaved flow 
predicted in this situtation. 

Predicted velocity contours at the exit of 
the diffuser are presented in Fig. 19. The flow is 
well behaved in this case with the velocity 
uniformly varying across the tunnel. Regions of 
flow separation or backflow are shown in each of 
the four corners. Figure 20 displays a predicted 
vertical velocity profile at the diffuser exit which 
also demonstrates the uniformly varying velocity 
across the tunnel here. 

The last case to be discussed is similar to 
the previous case in which the grid includes the 
contraction geometry, but in this case the 
turbulence model is activated. Figure 21 shows 
predicted velocity contours at a station halfway 
down the test section corresponding to the first 
station in the first case. As would be anticipated, 
the boundary layer is slightly thicker than in the 
first case, but the flow pattern is regular and 


Fig. 19 Predicted Vel. Contours, Diffuser Exit, 
Extended Grid, Laminar 


Fig. 20 Predicted Vel. Profile, Diffuser Exit, 
Extended Grid, Laminar 




Fig. 21 Predicted Vel. Contours, Test Section, 
Extended Grid, Turbulent 

symmetric around the tunnel walls. 

Figure 22 presents predicted velocity 
contours at the diffuser entrance. These contours 
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Fig. 22 Predicted Vel. Contour, Diffuser Entrance 
Extended Grid, Turbulent 


are also regular and symmetric and, in fact, are 
very similar to the contours at the station halfway 
down the test section which is 10 feet upstream. A 
vertical velocity profile is shown in Fig. 23 which 
also indicates a well behaved predicted flow 
pattern. 


Fig. 23 Predicted Vel. Profile, Diffuser Entrance 
Extended Grid, Turbulent 



Fig. 24 Predicted Vel. Profile, Diffuser Exit 
Extended Grid, Turbulent 


Figures 24 and 25 show predicted velocity 
contours and the vertical velocity profile, 
respectively, at the diffuser exit. The flow pattern 
is similar to the previous, laminar flow case, 
however, less separation or backflow is predicted in 
each of the four corners of the cross section. 


Fig. 25 Predicted Vel. Profile, Diffuser Exit 
Extended Grid, Turbulent 


In contrasting the measured velocity 
profiles with those from the computed results, 
several things are evident. First, the predicted 
boundary layer thickness for the case with the grid 
which does not include the contraction section is 
somewhat thin at the diffuser entrance being about 
6% of the tunnel height while both the extended 
grids agree better, at 9% for the laminar and 11% 
for the turbulent case. The measured boundary 
layer thickness is approximately 10% of the tunnel 
height. Furthermore, in contrasting the shapes of 
the velocity profiles at the diffuser exit, it is 
evident that the extended grid, laminar flow case 
shown in Fig. 20 most closely resembles the shape 
of the measured profile shown in Fig. 4, than do 
either of the two turbulent cases shown in Figs. 15 
and 25. 

CONCLUDING REMARKS 

The fundamental purpose of this paper was 
to numerically simulate the flow in the diffuser 
section of the IRT. It was found that the PARC3D 
computer code did a more acceptable job of 
predicting the flow details than did the previously 
used PEPSIG code. However, complete agreement 
with measured velocity values is yet to be realized. 
Several steps may be taken to improve this 
agreement: 

1. Further refinement of the grid structure is 
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3. Downstream turning vane effects need to be 
incorporated into the simulation. These will help to 
generate the asymmetry conditions seen in the 
measured flowfield values but not present in the 
predictions. 

4. Flow angularity measurements need to be made 
in the contraction section and/or the test section to 
verify that the uniform inflow boundary condition 
is valid. 

Besides the above, it is recommended that 
the numerical simulation of the IRT include the 
following: 

1. Numerical solution of the energy equation so 
that tunnel temperature distributins could also be 
predicted and contrasted to measured values. 

2. A particle trajectory code should be employed to 
predict the pattern of the icing cloud spray. 
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warranted. In particular, the wall layers require 
more grid points than were used in the current 
study. 

2. An investigation needs to be performed to 
ascertain that the predicted results are nongrid 
dependent. 
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